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ABSTRACT 


The present study was conducted to perform a molecular comparison of Pulex irritans based on the mitochondrial 
genome in four different climatic regions including the Caspian Sea region, a mountainous region, Persian Gulf region, 
and the Central Desert region, and based on nuclear ribosome genome in the west and northwestern Iran. A total of 
1937 adult flea samples were collected including 1019 P. irritans (52.61%) and 918 Ctenocephalides canis (47.39%) from 
various hosts including humans (14.1%), sheep (22%), goats (33.5%), dogs (25.6%) and houses (6.7%) between April 
2018 and May 2019. The samples collected from different hosts had similar morphological characteristics. However, 
there were slight differences based on mitochondrial markers and nuclear ribosomal markers in the study populations. 
The results from the phylogenetic tree based on three nuclear ribosome and mitochondrial markers showed that despite 
the slight differences in this sequence of different hosts and cities, all samples from different regions are in the same 
phylogeny. The results of ribosomal and mitochondrial genome analysis showed that these pieces are useful for demon- 


strating intraspecific similarity, and differentiation at species level and genus of P. irritans. 
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Introduction 


pee irritans belongs to Siphonaptera as a rela- 
tively-small order of wingless holometabolous 
insects, known as ectoparasites in their adult stage [1]. 

Certain flea species can transmit pathogens such 
as Yersinia pestis, Rickettsia typhi, and Bartonella 
henselae through their bites, feces and saliva. Biting 
can also cause severe itching and skin infections [2, 3, 
4]. P. irritans can cause serious health complications 
and transmit diseases in many parts of the world [5]. 
Global warming, international trade, travel, and pop- 
ulation growth influence the epidemiology of these 
infections [6, 7, 8]. 

Cladistic analysis using morphological features 
including head, chest, and abdomen have been con- 
ducted to classify various fleas, nevertheless, accurate- 
ly identifying certain flea species is difficult owing to 
variations in their morphological characteristics. Mo- 
lecular markers are therefore used to accurately iden- 
tify different flea species and investigate their phyloge- 
netic relationships [9, 10, 11, 12]. 

Molecular markers can be used to effectively char- 
acterize and specify an insect owing to their high sta- 
bility and conserved nature; for instance, ribosomal 
DNA (rDNA) with its subunits such as 28s, 18s, 5.8s, 
ITS1 and ITS2 is commonly used to identify insect 
species [9, 10, 11]. As the largest protein-encoding 
gene in the mitochondria of metazoan organisms, the 
cytochrome oxidase subunit I (COXI) gene has been 
extensively used to perform phylogenetic research 
and identify species and their differences [13, 14]. 

Accurately identifying flea species is crucial given 
that they constitute a vector of dangerous pathogens 
and thus a serious threat to public health. This study 
employed ITS1, ITS2, and COX1 as molecular mark- 
ers to investigate the molecular characteristics of dif- 
ferent populations of P. irritans and their phylogenetic 
relationships in four areas in Iran. 


Results 


After morphological examinations using reliable 
identification keys for Iranian fleas, 1019 (52.61%) out 
of 1937 samples were identified as P irritans (Table 1). 
Furthermore, the amplicons obtained from PCR for 
ITS1, ITS2 (rDNA), and COX1 (mtDNA) were ap- 
proximately 1000 bp, 500 bp and 700 bp, respectively. 
The Clustal Omega based comparison of the ITS1 se- 
quence in P irritans showed a similarity of 99.59% and 
a molecular diversity of 0.41% between the sequenc- 
es in ten different regions (Figure 1). A 100% simi- 
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larity was observed among the sequences of Sanan- 
daj (MN684797.1), Khorramabad (MN684787.1), 
Urmia (MN684784.1), Bahar (MN684777.1), and 
Kermanshah(MN 684776.1). These sequences 
were also similar (99.79%) to those of Hamadan 
(MN684779.1), Kuhdasht (MN684790.1), Kamyaran 
(MN684786.1), Gilangharb(MN684775.1), and Ma- 
habad (MN684778.1) with four nucleotide differenc- 
es. Four nucleotide polymorphisms were observed at 
positions 152, 165, 193 and 279 within the ITS1 spac- 
er in the ten populations. Sequencing the ITS1 frag- 
ment showed three 99-bp replicas in the first half of 
the sequence beginning at positions 146, 246, and 352 
and ending at position 431. Comparing the ITS2 gene 
sequences in the ten regions showed a 100% similar- 
ity between the sequences (Figure 2). Comparing the 
COX1 sequence of the nuclear genome (mtDNA) in 
the four different geographical zones (accession ID: 
MN173748-761) found all the species to belong to P 
irritans and showed a 99.86% similarity between the 
sequences (Figure 3). 

Multiple alignments of amino acids in COX1 
showed a 99.86% similarity and only one difference at 
position 54 of the methionine isolate and other isole- 
ucine isolates in Hamadan (Figure 4). Two subclades 
of the phylogenetic tree based on the similarities be- 
tween COX1 sequences in the present study and reg- 
istered sequences in GenBank comprised subclade 
Al, including P. irritans of the present and previous 
studies conducted in Iran and that isolated in China, 
Turkey, Jerusalem, Spain and Croatia. Subclade A2 
included Ctenocephalides felis in Iran and C. canis in 
Turkey. Subclades A3 and A4 included C. orientis and 
C. canis, respectively. Iranian subclade B also included 
Nosopsyllus fasciatus and Xenopsylla cheopis (Figure 
5). 

Three subclades of the phylogenetic tree based 
on the similarities between ITS1 sequences in the 
present study and registered sequences in GenBank 
comprised subclade A, including the sequences of 
the present and previous studies conducted in Iran. 
Subclades B and C included C. canis and N. fasciatus, 
respectively (Figure 6). 

Four subclades of the phylogenetic tree based 
on the similarities between ITS2 sequences in the 
present study and registered sequences in GenBank 
comprised subclade C, including the sequences of the 
present study. Subclades A and B included C. canis 
and X. cheopis, respectively . Subclade D also included 
N. fasciatus (Figure 7). 

The morphological characteristics of P. irritans 
were consistent with the molecular findings. Analyz- 
ing the partial COX1 gene can help identify P. irri- 
tans species and assess their intraspecific similarity. 
ITS1 and ITS2 sequences could also be used as useful 
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Figure 1. 
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Alignments of partial ribosomal DNA (ITS1) sequences in Pulex irritans isolated from sheep, human, goats, and dogs in the west and 
northwest of Iran (Kho(sheep), Bah(sheep), Ker(human), San(goat), Urm(sheep), Kuh(human), Kam(human), Ham(sheep), Ghi(goat), 
and Mah(dog) with accession nos. MN684787.1, MN684777.1, MN 684776.1, MN684797.1, MN684784.1, MN684790.1, MN684786.1, 
MN684779.1, MN684775.1, MN684778.1) Nucleotides in horizontal boxes display repeated units, and vertical boxes indicates the 
polymorphic site. Bah: Bahar; Ghi: Gilan-e Gharb; Ham: Hamedan; Kam: Kamiyaran; Ker: Kermanshah; Kho:Khorramabad; Kuh: 


Kuhdasht; Mah: Mahabad; San:Sanandaj; Urm: Urmia. 


markers for species-level differentiation and intra- 
specific similarity determination. The COX1 gene is, 
however, more efficient than the ITS1 and ITS2 genes 
in detecting genus, species and intra-species similari- 
ties. The diversity at the genus level and even between 
members of the same genus caused by high levels of 
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replacement in the Internal transcribed spacer coun- 
teracted their beneficial effects. Given this diversity, 
the sequences of the heterogenic versions of ITS frag- 
ments are determined in cloning plasmids and then 
specifically, which requires laborious and costly stud- 
ies[15]. 
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Figure 2. 
Alignments of partial ribosomal DNA (ITS2) sequences in Pulex irritans isolated from sheep, human, goats, and dogs in the west 


and northwest of Iran (Mah(human), San(human), Ham(sheep), 


Bah(human), Ghi(goat), Kho(human), Kuh(sheep), Ker(dog), 


Kam(dog), and Urm(dog ) with accession nos. MN684803.1, MN684801.1, MN 684800.1, MN684796.1, MN684792.1, MN684794.1, 
MN684793.1, MN684792.1, MN684781.1, MN684780.1). Bah: Bahar; Ghi: Gilan-e Gharb; Ham: Hamedan; Kam: Kamiyaran; Ker: 
Kermanshah; Kho:Khorramabad; Kuh: Kuhdasht; San:Sanandaj; Urm: Urmia. 


Discussion 


Molecular methods have enabled the identi- 
fication of different ectoparasite species with high 
morphological similarities [16]. The present study 
employed phylogenetic and molecular approaches 
to compare P. irritans in four different geographical 
zones of Iran. The morphological characteristics of 
P. irritans were found to be consistent with those ob- 
tained in previous studies conducted in Iran to deter- 
mine the flea fauna in different hosts and cities. [17, 
4, 18]. 

The host specificity can affect the intraspecific 
genetic diversity given the high levels of intraspecific 
genetic variations in general parasite species, which 
cause their infestation of broad ranges of hosts [19]. 
Research suggests no specific hosts for fleas, as they 
prefer diverse types of the host, which demonstrates 
intraspecific genetic diversity. Hornok et al. (2018) 
reported diversity in synanthropic flea species such 
as C. canis and P. irritans by investigating mitochon- 
drial sequences [20]. The present morphological data 
showed no differences between the P. irritans speci- 
mens collected from different geographical zones. The 
results of the morphological investigation were also 
consistent with those of the molecular exploration. 
The present study reported a 100% similarity between 
the nucleotide sequence of COX 1 gene in the P ir- 
ritans and the sequence of this gene in the P. irritans 
in similar studies conducted in Iran. [MF380389.1, 
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MF380390.1, MF380391.1 (NCBI)]. Similarities 
to the sequences in Spain (LT797470.1) and Chi- 
na (MF000666.1) were also reported as 99.85% and 
99.55%, respectively. Furthermore, the present find- 
ings were consistent with those obtained by Seyyedza- 
deh et al. (2018), who reported no differences between 
the isolates of C. canis in different areas of West Azer- 
baijan, Iran [21] and with those obtained by Zurita et 
al. (2019), who observed no significant differences in 
morphological data between P. irritans in Spain and 
Argentina. In addition, Zurita et al. (2019) reported 
a significant intraspecific similarity between the two 
populations based on mitochondrial genes [22]. Hor- 
nok et al. (2018) observed no morphological differ- 
ences between human and wild carnivorous P. irritans 
specimens in Hungary and Croatia [20], which is con- 
sistent with the present results. In contrast, Krasnov 
et al. (2015) found morphological differences between 
the flea species isolated from different hosts in differ- 
ent geographic areas to suggest high levels of genetic 
diversity [23]. 

Evidence suggests ITS1 and ITS2 constitute ap- 
propriate molecular markers for analyzing phyloge- 
netic relationships in fleas at their species level [24, 
11]. The present research observed a 99.59% similari- 
ty between the nucleotide sequences of the ITS1 spac- 
er, in five provinces in Iran. Analyzing the nucleotide 
sequence of the ITS1 spacer, also showed three 99-bp 
replicas in its first half, which is consistent with the 
results of Ghawami et al. (2018), who compared the 
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Alignments of partial mitochondrial DNA (COX1) sequences of Pulex irritans among populations from this study (Urm(human), 
San(sheep), Sar(sheep), Mah(goat), Kuh(sheep), Kho(human), Esf(human), Kam(sheep), Ghi(sheep), Ahv(goat), Ker(human), Kerm 
(sheep), Ham (human), and Bah(human) with accession nos. MN173761.1, MN173760.1, MN 173759, MN173758.1, MN173757.1, 
MN173756.1, MN173755.1, MN173754.1, MN173752.1, MN173750.1, MN173749.1, MN173748.1, MN173753.1, MN173751.1 ). The 
vertical boxes indicate the polymorphic sites. Bah: Bahar; Ghi: Gilan-e Gharb; Ham: Hamedan; Kam: Kamiyaran; Ker: Kermanshah; 


Kho:Khorramabad; Kuh: Kuhdasht; San:Sanandaj; Urm: Urmia. 


nucleotide sequence of the ITS1 spacer, in P irritans 
in two geographic areas and reported only one nucle- 
otide difference with a 99.85% similarity [25]. In line 
with the present results, those obtained by Vobis et al. 


(2004) suggested a relatively-constant nucleotide se- 
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quence of the ITS1 spacer, in different C. felis popu- 
lations [11]. 

A 100% similarity was observed between the ITS2 
nucleotide sequence of different samples. Zurita et al. 
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Figure 4. 
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Sequence alignment of COX1 amino acids for Pulex irritans in different geographical areas of Iran (Kerm (sheep), Bah (human), 
Ham (human), Ker(human), Ahv(goat), Ghi(sheep), Kam(sheep), Esf(human), Kho(human), Kuh(sheep), Mah(goat), Sar(sheep), 
San(sheep), Urm(human), with accession nos MN173748.1, MN173751.1, MN173753.1, MN173749.1, MN173750.1, MN173752.1, 
MN173754.1, MN173755.1, MN173756.1, MN173757.1, MN173758.1, MN 173759.1, MN173760.1, MN173761.1). The vertical boxes 
indicate the polymorphic sites. Ahv: Ahvaz; Bah: Bahar; ESF: Esfarayen; Ghi: Gilan-e Gharb; Ham: Hamedan; Kam: Kamiyaran; Ker: 
Kermanshah; Kerm: Kerman; Kho:Khorramabad; Kuh: Kuhdasht; San:Sanandaj; Urm: Urmia. 


(2019) reported an intraspecific variation of 90.1%- 
100% in four C. canis clones based on the ITS2 spac- 
er,, which is consistent with the present results. Vobis 
et al. (2004) also used the ITS2 sequence to identify 
different flea species and determine their intra-species 
differences [11]. P. irritans and C. canis belonging to 
the Pulicidae family constitute the main flea genera 
infesting human and cattle in the west and north- 
west of Iran. No significant differences were observed 
among the samples collected from different hosts. In- 
significant differences were found between the study 
populations in terms of mitochondrial markers and 
the nuclear ribosome. The results of plotting the phy- 
logenetic tree showed all the samples collected from 
different regions in the same branch despite their neg- 
ligible sequence differences. Analyzing mitochondrial 
genome also showed that this sequence is more useful 
than the nuclear ribosomal genome in illustrating in- 
tra-species similarity and differentiation at genus and 
species levels. 


Materials & Methods 


Location 


The Iranian plateau is located in the northern hemisphere 
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at 32.4279 °N, 53.6880 °E. The four study geographical zones in 
Iran included region 1: Caspian Sea (temperature: 8.00-26.00 
°C, annual precipitation: 400-1500 mm), region 2: Mountainous 
area (temperature: -5.00-29.00 °C, annual precipitation: 200-500 
mm), region 3: Persian Gulf (temperature: 12.60-35.00 °C, annu- 
al precipitation: 200-300 mm) and region 4: The Central Desert 
(temperature: -4.00-44.00 °C, annual rainfall: below 100 mm). Ac- 
cording to meteorological data, Iran is climatically divided into 
three zones, i.e. warm-dry, cold, temperate-humid, and warm-hu- 
mid [26]. Four different geographical zones in this country from 
which adult fleas were collected included Kermanshah, Kordes- 
tan, Hamadan, Lorestan, West Azerbaijan, Khuzestan, Khorasan, 
Mazandaran, and Kerman (Figure 8). 


Sampling 

A total of 1937 flea samples were directly isolated from the 
host body and houses using light traps and human-baited traps. 
They were transferred in 70% ethanol to the parasitology laborato- 
ry of the Faculty of Veterinary Medicine. They were then cleaned 
with distilled water several times, immersed in 5% potassium hy- 
droxide for 24 hours, re-washed with distilled water, and dehy- 
drated with a graded series of alcohol (50%-96%(. The fleas were 
ultimately identified by observing the slides with an optical micro- 
scope according to diagnostic keys [27, 28]. 


Morphological identification 


All the P irritans samples collected from different provinces 
and hosts were morphologically characterized by the following 
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Figure 5. 

The phylogenetic tree of Pulex irritans based on COX1. The molecular phylogenetic analysis was performed and the evolutionary his- 
tory was inferred by using the Maximum Likelihood method. The log likelihood (-1898.1614) of the phylogenetic tree was the highest. 
The percentage of the tree including a cluster of the relevant taxa is shown beside the branches whose length was measured as the num- 
ber of substitutions per site on this scale drawing of the tree. The analysis involved 44 nucleotide sequences. All the positions containing 
gaps and missing data were eliminated. The final dataset contained 448 positions. Evolutionary analyses were conducted in MEGA6. 
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Figure 6. 

The phylogenetic tree of Pulex irritans based on ITS1. The molecular phylogenetic analysis was performed and the evolutionary 
history obtained using maximum likelihood. The log likelihood (-1703.2418) of the phylogenetic tree was the highest. The percentage 
of the tree including a cluster of the relevant taxa is shown beside the branches whose length was measured as the number of substi- 
tutions per site on this scale drawing of the tree. The analysis involved 19 nucleotide sequences. All the positions containing gaps and 
missing data were eliminated. The final dataset contained 436 positions. Evolutionary analyses were conducted in MEGA6. 
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Figure 7. 

The phylogenetic tree of Pulex irritans based on ITS2. The molecular phylogenetic analysis was performed and the evolutionary history 
obtained using maximum likelihood. The log likelihood (-1272.22468) of the phylogenetic tree was the highest. The percentage of the 
tree including a cluster of the relevant taxa is shown beside the branches whose length was measured as the number of substitutions per 
site on this scale drawing of the tree. The analysis involved 26 nucleotide sequences. All the positions containing gaps and missing data 
were eliminated. The final dataset contained 291 positions. Evolutionary analyses were conducted in MEGA6. 
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Figure 8. 
Provinces of the study in the four geo- 
graphical areas included Kermanshah, 
Kordestan, Hamadan, Lorestan, West 
Azerbaijan, Khuzestan, Khorasan, 
Mazandaran and Kerman. 


features: the round-headed anterior margin with a pair of eyes 
and no protrusion, absence of genal and pronotal ctenidia, asym- 
metrical club of antenna, strong setae under the eyes, row of small 
spines mass on the inner surface of the posterior side of the coxa 
in females 7-10 and in males 8-12, and no pleural rod on the me- 
sopleural segment [27]. Comparing the samples collected from 
different provinces and hosts based on morphological character- 
istics identified 52.61% to belong to P irritans species despite their 
negligible differences. 


Molecular analysis 


The samples were individually crushed with a sterile scalpel, 
poured into a sterile microtube and their total DNA was extract- 
ed using DNA extraction kits (MBST, Tehran, Iran) according 
to the manufacturer’s instructions. The extracted DNA samples 
were stored in sterile microtubes at -20 °C until use. DNA was 
extracted from the samples isolated from different hosts in dif- 
ferent areas. The specific primers used by Vobis et al (2004). were 
employed to identify P irritans and 1000-bp ITS1 and 500-bp 
ITS2 fragments in the rDNA of P irritans. Primers used by Fol- 
mer et al. (1994) were used to amplify the mitochondrial COX1 
gene fleas P. irritans (Table 2). The result of this proliferation was 
observed as a 700 bp band [29]. To determine molecular differ- 
ences between the species isolated from different hosts and areas, 
PCR was performed in a total volume of 50 ul, containing a 5-ul 
DNA template, a5-ul 10x PCR buffer, 1-ul Taq Polymerase, 1 ul of 
each primer (20 uM), 1 ul of dNTP (100 uM) and 4 ul of MgCl2 
(50 mM). All these materials were procured from SinaClon, Iran. 
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The thermal cycle of the PCR on fragments ITS1 and ITS2 was 
as follows: incubation at 94 °C for 5 minutes to denature dou- 
ble-stranded DNA followed by 30 cycles of one-minute denatur- 
ation cycles at 94 °C, annealing at 54 °C for 1 minute, extension 
at 72 °C for 1 minute and an additional extension at 72 °C for 5 
minutes. The COX1 gene also underwent replication at 95 °C for 
5 minutes followed by 35 cycles of denaturation cycles at 95 °C for 
1 minute, annealing at 55 °C for 45 seconds, extension at 72 °C 
for 45 seconds, and final extension at 72 °C for 10 minutes. The 
PCR products were electrophoresed on 1.5% agarose gel, stained 
with a safe stain, visualized under UV light and then purified for 
sequencing. After purifying the PCR products, 34 samples, in- 
cluding 10 positive samples for ITS1, 10 positive for ITS2, and 14 
positive for COX1 sequence, collected from different hosts were 
transferred to Takapouzist Company for sequencing. 


Data analysis 


A total of 34 PCR products were transferred to Takapouzist 
for sequencing. An accession number was obtained by online 
recording the sequencing results on NCBI by the host and geo- 
graphical area. The data were analyzed using blasting sequences 
on NCBI and plotting the phylogenetic tree in Mega 6 using max- 
imum likelihood and bootstrapping (1000 replicates). The cluster- 
ing method proposed by Zurita et al.(2019) was also conducted. 
The required sequences were ultimately aligned and compared 
using EMBOSS Needle and Clustal Omega to evaluate the sim- 
ilarity percentage of the nucleotide sequences of mitochondrial 
(mtDNA) and ribosomal genomes (rDNA) in different provinces 
and hosts. 


Table 2. 
Nucleotide sequence and specificity of the primers used. 
Annealing Fragment length 
Gene Primer nucleotide sequence temperature a s 1 Reference 
(°C) (base pair; bp) 
For: GTA CAC ACC GCC CGT GCG TAC T A 
a Rev: GCT GCG TTC TTC ATC GAC CC as ae [23] 
For: GGG TCG ATG AAG AAC GCA GC ò 
1152 Rev: GCG CAC ATG CTA GAC TCC GTGGTT CAA G al 200bp [29] 
For: GGT CAA CAA ATC ATA AAGATA TTG G 
COX1 a 55°C 700 bp [30] 


Rev: GAA GGG TCA AAG AAT GAT GT 
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